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ABSTRACT 

An analytical model is developed for the mass function of cold dark matter subhalos at the time 
of accretion and for the distribution of their accretion times. Our model is based on the model of 
Zhao et al. (2009) for the median assembly histories of dark matter halos, combined with a simple 
log-normal distribution to describe the scatter in the main-branch mass at a given time for halos of 
the same final mass. Our model is simple, and can be used to predict the un-evolved subhalo mass 
function, the mass function of subhalos accreted at a given time, the accretion-time distribution of 
subhalos of a given initial mass, and the frequency of major mergers as a function of time. We test our 
model using high-resolution cosmological Y-body simulations, and find that our model predictions 
match the simulation results remarkably well. Finally, we discuss the implications of our model for the 
evolution of subhalos in their hosts and for the construction of a self-consistent model to link galaxies 
and dark matter halos at different cosmic times. 

Subject headings: cosmology: dark matter halos - galaxies: formation - galaxies: halos 



1. INTRODUCTION 

In the current Cold Dark Matter (hereafter CDM) 
paradigm of structure formation, a key concept in the 
build-up of structure in the universe is the hierar- 
chical formation of dark matter halos. Galaxies and 
other luminous objects are assumed to form by cool- 
ing and condensation of baryons w ithin these halos (see 
lMo~ van den Bosch fc White 1 120101 for an overview). In 
this scenario, a detailed understanding of the formation 
and structure of dark matter halos is of fundamental im- 
portance for predicting the properties of luminous ob- 
jects, such as galaxies and clusters of galaxies. 

The formation history of a CDM halo is conve- 
niently represented by its merger tree, which de- 
scribes how its progenitors merge and accrete during 
its entire formation history. For a given cosmologi- 
cal model, such merger trees can be constructed ei- 
ther from N-body simulations or from Monte-Carlo re- 
alizations based on the extended Press-Schechter (PS 
formalism (IPress fe Schechterl U974t IBond et all 119911 



BowerillfMlLacev fc Cole 111993b ISheth. Mo fc Tormanl 
20011) . In recent years, much effort has been made 
to characterize and understand the statistical proper- 
ties of halo formation in a CDM cosmogony. One 
particular aspect of the halo formation process is the 
existence of a subhalo population, which is produced 
by the accretio n and survival of p rogenitors at vari- 
ous times (e.g . iKlvpin et all 119991: iMoore et al.l 119991: 



KravtS0^ i e^_alj_j£UU2Lj^aye^iij l^uu^J^^iua^ai. 

20041: Ivan den Bosch et al.l 120051: iDiemand et al.l 12007 
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[20091: lAngulo et al.l [20091: iLi et all [2001 ). Since galaxies 
may form at the centers of these progenitor s and merge 
into t he final halo along with their hosts (e.g. lKang et all 
2005), the statistical properties of the subhalo popula- 
tion are expected to be closely linked to those of satellite 
galaxies. One of the basic properties of the subhalo pop- 
ulation is the mass function of the progenitors of subha- 
los (i.e., the masses of the subhalos at their mo ment of 
accretion). Following Ivan den Bosch etall (pOOl . we re- 
fer to this mass function as the un-evolved subhalo mass 
function, to distinguish it from the evolved subhalo mass 
function that refers to the present day masses of dark 
matter subhalos (see f}3] for details). A number of re- 
cent investigations have used these subhalo mass func- 
tions in an attempt to characterize the galaxy-dark mat- 
ter connection across cosmic ti mes (e.g.lVale fc Os triker 
200 l [2001 iZheng et al.ll2005l [2007t iConroy et al l | 2006l 
2007t IConrov fc Wechslerl 120091: lYang et all 120091 poTl) 
Li et " all 12009) iMoster et all 120101: iBehroozi et aUboiot 
Wang fc Jingfeo lOt IWetzel fc WhitefeOlOHNeistein et all 



Wang & Jmg ^UlU; Wetzel & White 
201lHAvila-Reese fc Firmanill201lD . 



So far the subhalo mass function has been studied 
with Y-body simulations and Monte- Carlo realizations 
of th e extended PS forma li sm (e .g. ISheth fc Lemsor 
1991 



ISomerville fc Kolatti [T9991 ICole et all 
van den Bosch et al.l 120051 iGiocoli et al ' 



Cole et al.ll200a IParkinson et al.1120081: iFakhouri fc Mai 
20081: IFakhouri et al.l 120101 ) . in this paper we show 



that a simple analytical model can be constructed to 
describe not only the mass distribution of subhalos but 
also the distribution of their accretion times. We use 
Y-body simulations to demonstrate that the model is 
remarkably accurate. The model is not only simple to 
implement, but also provides important insights into 
the formation and evolution of dark matter halos and 
subhalos. Furthermore, as we briefly discuss in ^ our 
model also provides a self-consistent way to link galaxies 
and dark matter halos at different redshifts. 

The structure of the paper is organized as follows. In 
Section [5] we outline the simulations based on which we 
will test our model. In Section [3] we describe our model. 
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In Section 2] we present our model predictions for the 
conditional mass function of subhalos, the major merger 
rates, and the subhalo mass function. These model pre- 
dictions are tested against TV-body simulation results in 
Section|4l Finally, in Section[5]we discuss the universality 
of our model for other models of structure formation, and 
discus how our model can be used to study the evolution 
of subhalos in their hosts and to construct self-consistent 
models that link galaxies and dark matter halos across 
cosmic times. 

Throughout this paper, we use 'In' to denote natural 
logarithm and 'log' to denote the 10-based logarithm. 

2. THE SIMULATIONS 

Before presenting our model, let us first describe briefly 
the the iV-body simulations to be used to test the model. 

We use two different A-body simulations that assume 
the same cosmology but use different box sizes (mass 
resolutions). Both simulations were car ried out using 
the m a ssively parallel GADGET2 code (Springcl ct al. 
I2001al [2005h . The simulations evolved 1024 3 dark 
matter particles in periodic boxes of 100 /i~ 4 Mpc and 
300 /i _1 Mpc on a side, respectively, from redshift z = 
100 to the present epoch (z = 0). The particle 
masses and softening lengths are, respectively, 6.93 x 
10 7 h- l M e and 2.25 /i^kpc for the 100 /i^Mpc box 
simulation, and 1.87 x 10 9 h~ 1 M Q and 6.75 /i _1 kpc for 
the 300 h~ 1 Mpc box simulation. The cosmological pa- 
rameters usedTn_rtie_^imulations are based on those pub- 
lished in ElmkiixjEII] ((111): fl m = ft dm + ft b = 0.258, 
fib = 0.044, h = 0.719, fl A = 0.742, n = 0.963, and 
o 3 = 0.796. Unless stated otherwise, our model predic- 
tions are also for the ACDM cosmology with this par- 
ticular set of parameters. For both simulations, a total 
of 100 outputs in equal log(l + z) interval were made, 
starting from z — 50 (for 300 /i _1 Mpc box) and z = 20 
(for 100 h^Mpc box) to z = 0. 

Dark matter halos were identified from the simula- 
tions at each outp ut using the stand ard friends-of-friends 
(FOF) algorithm (jDavis et al.ll985f ) with a linking length 
of 0.2 times the mean interparticle separation. Here we 
keep all halos with at least 20 particles. Based on halos 
at dif ferent outputs, halo merger trees were constructed 
(see lLacev fc"CoIe1ll993l) . A halo in an earlier output is 
considered to be a progenitor of the present halo if more 
than half of its particles are found in the present halo. 
The main branch of a merger tree is defined to consist of 
all the progenitors one goes through as one climbs from 
the bottom to the top, choosing always the most mas- 
sive branch at every branching point. These progenitors 
are referred to as the main-branch progenitors, and the 
time dependence of the main branch mass is referred to 
as the assembly history. In our analysis based on the 
300 /i _1 Mpc box simulation, we randomly choose about 
200 trees, from the total merger tree catalogue, to sam- 
ple each of the two massive bins at Mh ~ 10 145 /t _1 MQ 
and M h ~ 10 14 /i _1 M Q , and about 2000 trees to sam- 
ple each of the low-mass bins at Mh ~ 10 13 ft _1 M0 
and M h ~ 10 120 fr -1 M , with bin widths AlogM^ ~ 
0.7, 0.1, 0.1, 0.02, respectively. Our tests later are mainly 
based on the 300 /i~ 4 Mpc box simulation. However, 
in many cases we also use halo merger trees obtained 
from the 100 /i _1 Mpc box simulation to achieve bet- 



ter mass resolution. Specifically, about 7, 36, 400 and 
2000 trees are selected from this simulation to sample 
the merger histories for halos with Mh ~ 10 14 5 h~ 1 M & , 
~ 10 14 /i" 1 M , ~ 10 130 /i" 1 M Q and ~ 10 120 /i" 1 M , 
with bin widths A log Mh ~ 0.5, 0.5, 0.5, 0.3, respectively. 

3. THE MODEL 

Now we come back to our modeling of the accretion 
of subhalos. We want to obtain the distribution of dark 
matter subhalos with respect to their mass at accretion, 
m a , and their accretion redshift, z a , in a host halo of 
mass Mo at redshift zq. For convenience we use 

s a = a 2 a = a 2 (m a ); S = a 2 = <j 2 (M q ) (1) 

to label the masses, m a and Mq, and 



S a = S c (z a ); S = S c (z ) 



(2) 



to label the redshifts z a and Zq. Here a(M) is the vari- 
ance of the linear density field at z = on mass scale M, 
and 6 c (z) ps 1.686/D(z) [with D(z) the linear grow fac- 
tor normalized to unity at z — 0] is the critical density 
of spherical collapse at redshift z. We write the mean 
number of subhalos of mass m a accreted at redshift z a 
in a host halo (Mo, So) as 

d 2 N a =M a (s |5 ,(5o)dlnm a dln(l + z a ) . (3) 

The mean mass, M(z), of the main branch halos for 
all (Mq, So) ~ halos is in general a monotonically decreas- 
ing function of redshift (e.g.. lAvila- Reese et al.l 119981 : 
Ivan den Bosch et al.ll2002HWechsler et alJ^OO^ Tand we 
write M a = M(z a ). Hence, for given (M , So) we can use 
M a as a time-variable. Denote by F(s a , S a \So, So; M a ) 
the mean fraction of the total mass accreted in the 'time- 
interval' [M a — dM a ,M a ] that is in halos of mass m a . 
We can write 

d 2 N a = —F(s a ,S a \So,So;M a )dM a ds a . (4) 
m a 

It then follows that 

Kf t~I r in r TT >i &M a d In S a 

dm a din d a dln(l + z a ) 

Note that ds a /dm a is determined by the perturbation 
power spectrum, dlnc5 a /dln(l + z a ) by the linear growth 
factor, and dM(z a )/dhi5 a by the mean halo assembly 
history. 

We can integrate Af a over the mean mass assembly 
history to obtain the so-called un-evolved subhalo mass 
function: 



dN a 
din 777 



Af a (sa,S a \S 0l S )dln(l + z a ) 



Mn 





F(s a ,S a \S ,So;M a )-^dM a . (6) 
dm a 

This function describes the distribution of the masses at 
accretion of all subhalos accreted into the main-branch 
of the merger history of the (Mo, So) host halo. Thus 
for a given cosmology, one can obtain both Af a and 
the un-evolved subhalo mass function once models for 
J~(s a , S a \So, So', M a ) and for the mean halo assembly his- 
tory are adopted. In general T can be obtained using 
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halo merger trees constructed either from numerical sim- 
ulations or from analytical models, such as the extended 
PS formalism. Here we develop a simple analytical model 
based on the statistical properties of halo assembly his- 
tories, instead of on the full merger trees. 

Consider a halo whose main-branch mass is M a at z a . 
Since we are modeling the accretion of subhalos into the 
main branch of its merger tree, we must have that m a < 
m max witlfl 



T7?,, 



MIN(M o ,M /2). 



(7) 



Thus, a simple model for the mass fraction in (m a ,8 a )- 
progenitors to be accreted at z a may be written as 

$(s a ,5 a \S ,6 ;M a ) 

B^^^Fisa^SalSoJoiMa) if m a < m max 



otherwise , 



(8) 



where F(s a , 8 a \So, So; M a ) is the mass fraction in progen- 
itor halos of mass m a to be accreted at redshift z a . The 
normalization factor, 

/' °° 

B = / F(s \S ,5 ;M a )ds a , (9) 

JS(m m „) 

is the total mass fraction of all progenitors that can be 
accreted into the main branch. Now suppose that the 
distribution of M at z a is given by P(M a \So, So) then 



7"(s a) 5 a |S ,(5o;M Q ) 

$(sa, S a \S , S Q ; M a )P(M a \S ,S )dM a . 



(10) 



3.1. Models for F(s a ,8 a \S , So; M a ) 

Model I: According to the extended PS formalism, 
the fraction of mass of halo (Mi , zi) that is in progenitor 
halos of mass M2 at redshift Z2 > z\ can be written as 



f(S2,62\SM = 



1 S 2 -5 X 



2^ (S 2 - Sif/- 



■ exp 



(&2-&1Y 



2(S 2 



Si) 
(11 



(see lLacev fc Colelll993D . Thus, the simplest model is to 
assume F(s a , S a \S , S ; M a ) = f(s a ,5 a \S ,S ). However, 
f(s a , S a \So, So) is the mean for all (5*o, 5o)-halos, and so 
such an assumption completely ignores the scatter in the 
assembly history, i.e. the dependence on M a . A better 
approximation is to assume that F(s a , S a \So, So; M a ) is 
independent of So and So, i.e. the accretion properties at 
z a is determined entirely by the mass of the main-branch 
halo at that redshift regardless of where the halo will end 
up at z — 0. In this case, we may write 

\So, So; M a ) din 5 a 
= f[sa,5 a + dS a \S(M a ),5 a ] 
1 dS a 



^ [s a - 5(M„)]3/a ' 



(12) 



where the extended PS formula (ITT1) is used in the second 
equation. In what follows we will refer this model as 

5 Note that m a reflects the mass of the progenitor halo which, 
after accretion into the main-branch, increases the main-branch 
mass to M a 



Model I. Note that as we will illustrate in section B~Tl this 
model does not match well with the simulation results. 

Model II: Numerical simulations have shown that 
the PS formula (fTTj) is not accurate. Attempts have 
been made to come up wi t h better approximations 
(e.g. ISheth fc Tormenl l200l INeistein fc Dekel I [200l 
ICole et al.l 120081: IParkinson et al.l 12008ft . According to 
the em pirical modification proposed by IParkinson et al.1 
(2008]), we may write 

\So, So; M a ) dln<5 a 
1 d<5 a 



2^ [s a - 5(M a )]3/2 



G 



S a 



a(M a y a(M a ) 



(13) 



where G(x,y) = G a; 7l y 72 , with G = 0.57, 71 = 0.38 
and 72 = —0.01. This model will be referred to as Model 
II in what follows. 

Model III: Other than the above two models, we find 
by trial and error that the following simple modification 
of the extended PS formula provides a much more accu- 
rate model for F: 



F(sa,,S a \So,So;M a ) d\n5 a 
1 S a - S M 



2^ (s a - S M f/* 



exp 



2 1 



(<5q - 5 M ) 

2(s a - S M ) 



(14) 



where 8m corresponds to the redshift at which the main 
branch has the mass 

Mmax = MIN (.M a + TO max , M ) , (15) 

with M. a the median main-branch mass (not to be con- 
fused with the mean branch mass M a ) and 

Sm = o~ M = er 2 (Af max ) . 
This model will be referred to as Model III. 



(16) 



3.2. Halo Assembly History 

Following Z hao et all ()2009l ) , the median accretion rate 
of a (Mo, So) host halo at redshift z can be written as 

dlncr(M) 



din S c (z) 
Here 



5.85 



{w[S c (z), <t(M)] - p[5 c (z); S , a Q }} 



Scjz) -. n -dlncr(Af)/dlnM . 

a{M) 



and 

p[S c (z);So,o- ] 

= p(5 ;So,o-o) x Max 
with 

p(S a ;S ,ao) = 



(17) 
(18) 

(19) 



0,1 



log^z) - \ogS 



0.272/w{S ,a ) 
w(So,a a ) 



l + [w{8o,ao)/Af 



(20) 



One can obtain the median main-branch mass 
M a (z a \M , z ) by simply integrating Eq. (fTT)) from red- 
shift zq to z a - The solid curves in Fig. Q] correspond 
to the median assembly histories thus obtained for four 
different host halo masses, as indicated in each panel. 
For comparison, the thin, jagged curves correspond to 
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Fig. 1. — The thin solid curves in each panel are the assembly hist ories of 200 s imu lated halos of a given final mass as indicated in the 
panel. The solid cu rve is the predicted median assembly history by the Zhao ct al. (2009) model. The two dashed curves are the ±lcr-range 
given by equation (121 1 . 



mass accretion histories obtained from iV-body simula- 
tions. The results are shown for a selection of 200 main 
branch assembly histories for halos in each mass bin ex- 
tracted f rom the 300/i~ 1 Mp c box simulation. Clearly, the 
model of lZhao et al.l (|2009f) yields median mass assembly 
histories that are in excellent agreement with simulation 
results. 

In order to proceed we n eed t o make one important 
addition. The iZhao et al.l (|2009f ) model only gives the 
median assembly history. However, a complete model 
for the un-evolved subhalo distribution with respect to 
the mass at accretion and the redshift of accretion re- 
quires the full distribution function P(M a \So, So)- It is 
easy to understand that the dispersion in M a plays an 
important role. Since by definition the masses of subha- 
los are all smaller than that of the main progenitor, us- 
ing the mean assembly history would imply that at any 
given time no accreted subhalo can have a mass larger 
than the mean mass on the main branch. Clearly, when 
allowing for dispersion in the main-branch masses, this 



constraint is no longer present. The different panels of 
Fig. [5] show the distributions of the main-branch halo 
masses at different rcdshifts obtained from simulations 
(histograms). Results are shown for halos with redshift 
zero masses ~ 10 14 5 ft.~ 1 M Q (panels in the upper two 
rows) and ~ 10 13 /i _1 Mq (panels in the lower two rows), 
respectively. All the distributions are reasonably well de- 
scribed by a log-normal d istribution, with median given 
by the IZhao et all (|2009f ) model, and with a dispersion 
(in 10-bascd logarithm) given by 

a = 0.12 - 0.15 log(M a /M ) , (21) 

as shown by the solid curves in the figure. The lack 
of low-mass main-branch progenitors in the simulation 
apparent in the bottom right two panels is simply due 
to the mass limit in the iV-body simulation. For com- 
parison, the dashed lines in Fig. Q] show the ±lcr range 
given by Eq . (12T1) . together with the median given by the 
IZhao et ail ()2009l ) model (solid curves). 

An important advantage of the log-normal form for 
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Fig. 2. — The distribution of the main-branch mass at different redshifts obtained from ./V-body simulations (histograms). Results are 



shown for halos with final masses ~ 10 14 ' 5 h 'Mq (upp er two rows) and 
curves are the log- normal model described in Section 1.3.21 



10 1 



h Mq (lower two rows), respectively. The smooth 



the distribution P(M a \So, So) is that it is straightforward 
to compute the average accretion rate dM Q /dln(l + z a ) 
(see Eq. [5]). After all, for a log- normal distribution the 
average M a is related to the median Ai a according to 
M a = e( ln ( 1 °) cr ) I 2 M a - Thus, assu ming a log-norma l 
distribution, we can simply use the IZhao et all ([2009h 
model for the median assembly history to obtain the 
mean accretion rate. Fig. [3] shows the mass accretion 
rate d(M a /M )/dlog(l + z) for halos with different final 
masse s. The dashed cury e in each panel is the prediction 
of the IZhao et all (|2009t ) model of the median together 
with the log-normal distribution described above. 

In the simulations, especially for small halos at low red- 
shifts, there are time-steps over which the accretion rate 
is negative. This can come about due to, for instance, 
tidal stripping by neighboring structures, the loss of un- 
bound subhalos or unbound particles, the fragmentation 



of halo s, and the FOF-bridging problem. The Zha o et al.l 
(2009) model has taken such effects explicitly into ac- 
count via a correction factor in equation (j 19[) . In gen- 
eral, however, it is possible to have two different defi- 
nitions of the un-evolved subhalo population, one based 
on all subhalos that have entered the main branch at 
some point (but are not necessarily bound to or located 
within the final halo) , and the other based on those that 
have more than half of their particles ending up in the 
final halo. In this paper we adopt the first definition in 
which the accretion rate is determined by all subhalos 
that have been accreted onto the main branch at some 
point in time, regardless of whether they have left the 
main branch again or not. The open circles in Fig. [3] 
show the simulation results based on this definition. As 
one can see, the mean mass accretion rates predicted by 
the IZhao et ail (|2009T ) model match well the simulation 
results over a large range of redshift; the mismatch seen 
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Fig. 3. — The mean mass accretion rate of host halos with different final masses. The open circles in each panel show the rate in terms 
of dlog(l + z) obtained from TV-body simulations, normalized by the mass of the host ha lo, with error-bars obtained from 200 bootstrap 
resampling of the host halos. The da shed line is the model prediction of Zhao et al. (2009), while the dotted line shows the corrected mass 
accretion rate according to equation {22} . 



at high-z for low-mass halos is simply an artefact due to 
the limited numerical resolution of the simulation. The 
only significant discrepancy occurs for low-mass halos at 
low redshift where the predicted rate is lower than the 
simulated rat e. This dis c repan cy is mainly due to the 
fact that the iZhao et ail (|2009D model tries to account 
for the affects of tidal stripping and the presence of un- 
bound subhalos. For our definition of the un-evolved 
subhalo population, however, such correction is not re- 
quired when comp uting dM /dln(l+z). As discussed in 
IZhao et al.1 (poM ). this can be achieved by simply setting 

p[8 c (z);5 ,a ] =0 (22) 

in Eq. (fT"7|) . The dotted lines in Fig. [3] show the predic- 
tions based on Eq. (|2"2"|) . Clearly this simple modification 
works remarkably well, bringing the predictions in good 
agreement with the simulation results. 

If one is interested in using the second definition for 
the population of subhalos, one can still use our model 



but with Eq. replaced by Eq. (JT5J). In what follows, 
we present predictions of Model III based on the first 
definition, corresponding to Eq. (|2"2"j) . We have tested, 
though, that using the second definition instead yields 
results that are very similar, except when it concerns the 
subhalos accreted at low redshifts in low-mass host halos. 

4. TEST WITH iV-BODY SIMULATIONS 

In this section we use the -ZV-body simulations de- 
scribed in Section 2 to test our models for the distribution 
of subhalos with respect to their accretion redshift and 
their mass at accretion. 

4.1. The conditional mass function 

The first quantity we consider is the conditional mass 
function of subhalos in host halos of different masses 
(Eq. [TO]). From the halo merger trees constructed from 
the simulations, we measure the mass distribution of ac- 
creted subhalos within equal log(l + z)-bins centered at 
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rows). The symbols are results from N-body simulations, with error-bars obtained from 200 bootstrap resampling of the host halos. fn each 
panel, results are shown for both high- (100 /i -1 Mpc box; filled triangles) and low- (300 h~ 1 Mpc box; open circles) resolution simulations. 
In the high-resolution simulation, unbound particles are removed (see text for details). The long-dashed, dashed and solid lines show the 
predictions of Models I, II, and III respectively. Finally, the vertical dotted lines in each panel correspond to the mass limit of 100 particles 
in the high- (left line) and low- (right line) resolution simulations, respectively. Note that bin widths in different columns are different. 



0.42, 1.2, 2.4 and 4.2, respectively. The results are shown 
in Fig. [4] as open circles (for 300 /i _1 Mpc simulation 
box) and filled triangles (for 100ft. _1 Mpc simulation box) 
where the errorbars have been obtained from 200 boot- 
strap resamples of the population of host halos. Each col- 
umn shows results for host halos with a given mass, while 
each row shows the results at a given redshift. Note that, 
because of mass resolution and because of tidal effects, 
a fraction of the subhalos containing < 100 particles, es- 
pecially those near massive halos, are not gravitationally 
bound. This is the reason for the artificial upturn in 
the conditional mass function seen in the low-resolution 
simulation results at low redshift. In the literature this 
problem has been treated either by using only halos that 
contain large enough number of particles or by getting rid 
of unbound particles from halos. For instance, the SUB- 



FIND halo finder developed by iSpringel et all (|200ib[ ) 
tries to s eparate bound and u nbound structures in FOF 
halos. In lBenson et al.l (|2001l ). unbound particles are re- 
moved iteratively from each halo until the total energy of 
the halo becomes negative. In order to quantify the mag- 
nitude of this effect, we have also constructed halo mer ger 
trees using a method similar to iBenson et al.l (|200lD to 
calculate the number of bound particles per (sub)-halo. 
Halos with less than 20 bound particles are discarded 
from our sample. In Fig. |4l we show the results ob- 
tained from the 100 /i _1 Mpc simulation box in which we 
removed the halos with bound particles less than 20 (us- 
ing the method described above), together with the re- 
sults from the 300 /i _1 Mpc simulation box without such 
treatment. The dotted, vertical lines correspond to halos 
with 100 particles in each of the two simulations, and are 
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Fig. 5. — The major merger events per unit log(l + z) as a function of z. The symbols show results obtained from simulations, while the 
three curves in each panel show the predictions of Models I, II and III. For comparison, results for both the high-resolution (which extends 
to higher redshift; filled triangles) and low-resolution (open circles) simulations are shown for cases where the the statistic is sufficiently 
good. 



show for comparison. As one can see, the treatment of 
the unbound halos is able to remove the spurious upturn 
at the low-mass end. But for sure and to be conservative, 
the results presented below are all for subhalos contain- 
ing at least 100 particles. 

The long dashed lines in Fig. |4] show the predictions of 
Model I. This model does not match very well with the 
simulation results, especially for massive halos at low red- 
shifts. As shown in the upper-right panel of Fig. [4j for 
massive host halos this model significantly over-predicts 
the number of accreted low-mass subhalos and under- 
predicts that of accreted massive subhalos. Such discrep- 
ancy has already been noticed and extensively discussed 
in literature (e.^ 
iParkinson et al 



Sheth fc TormeiJ200alCole et al.H2008t 
2008> iNeistein et all boiOfl. An 



cm- 



irical modification was suggested by IParkinson et all 
2008), which is adopted in our Model II. The results 
obtained from this model arc shown in Fig. |4] as the 
dashed lines. As one can see, this modification suc- 



cessfully suppresses low-mass subhalos, but makes the 
under-prediction of massive subhalos worse. Finally, let 
us look at our Model III, the results of which are shown 
in Fig. H] as solid curves. Clearly, Model III matches the 
simulation results much better than either model I or II, 
especially for massive hosts. 

Note that all three models shown have adopted the cor- 
rection of Eq. (|22|) . If we use model III without this cor- 
rection, the model underpredicts the abundance of sub- 
halos (as defined using the first definition described in 
£13. 2p . especially for low- mass host halos at low redshifts. 
For completeness, we have also tested model III without 
the dispersion in halo assembly histories. As expected, 
not taking this dispersion into account yields conditional 
subhalo mass functions in poor agreement with the sim- 
ulation, especially at higher redshifts. Hence, we caution 
that it is important to properly account for the disper- 
sion in halo assembly histories. 
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4.2. Redshift distribution of major mergers 

With the conditional mass functions described above, 
it is straightforward to calculate the accretion rate of 
subhalos as a function of host halo mass and redshift. 
Here we focus on the rates of major mergers where the 
mass of the accreted subhalo is required to be m a > 
-M a /3. This mean rate in terms of redshift interval can 
be written as 

J^TT s = K(s a ,6 a \S ,5 )^. (23) 

dln(l + z Q ) J Ma /3 m a 

Note that since the major mergers are defined with re- 
spect to individual main branch masses M a , here we do 
not include the lognormal scatter (Eq. (|2"Tj) ) in perform- 
ing the integration. The predictions of Models I, II and 
III are shown in Fig. [5] as the long dashed, dashed and 
solid lines, respectively. To check the model predictions 
we calculate the same quantity directly from the N-hody 
simulations and the results are shown in Fig. [3] as sym- 
bols with error bars. For comparison, we show the results 
obtained from both the 300 /i _1 Mpc and 100 /i _1 Mpc 
simulation boxes. In the case of the more massive halos 
(upper two panels), however, we do not show the results 
from the 100/i _1 Mpc simulation box, since these are very 
noisy due to small number statistics. The 300 h~ 1 Mpc 
box, although having better statistics, has insufficient 
mass resolution to properly resolve the merger statistics 
of low mass halos at higher redshifts (as is clearly evident 
from the lower two panels) . Clearly, the results of model 
III are in very good agreement with the simulation re- 
sults, and fair much better than either Model I or Model 
II. In what follows, we focus only on the predictions of 
Model III. 

From the model for the major merger rates in halos of 
different masses and at different redshifts, we can also 
predict the characteristic redshift and the characteristic 
mass associated with the last major merger. By inte- 
grating Eq. (|23[) over redshift from to some z a , we can 
obtain the average number of major mergers, -/V ma j or , ex- 
pected in this redshift interval. We call the value of z a 



the characteristic redshift, and the main-branch mass at 
this redshift the characteristic mass of the last iV ma j or th 
major mergers. In Fig. [6] we show the characteristic red- 
shift (left panel) and characteristic mass (right panel) as 
functions of host-halo mass, for -/V ma j or = 1,2 and 3, cor- 
responding to the 1st, 2nd and 3rd last major merger, re- 
spectively. Since more massive halos assemble later (see 
Figd): their last major mergers, on average, occur at 
lower redshifts, as shown clea rly in the l eft-ha nd panel of 
Fig. [5] (see similar findings in lLi et all (|2007[ )). Interest- 
ingly, the characteristic mass of those last major mergers, 
in units of the host halo mass, is virtually independent 
of the host halo mass. This is a manifestation of the fact 
that all halos (in the mass range explored here) have 
average merger histories that are self-similar if they are 
allowed to be stretched or shortened along the time-axis. 
For example, it is interesting to notice that the average 
time interval between major mergers is roughly the same 
as the time during which the main branch mass increases 
by a constant factor, about 10°' 8 , quite independent of 
the host halo mass. 

In a recent paper based on the Millennium simula- 
tions, IFakhouri etail (|2010f ) have come up with a fitting 
formula which describes the subhalo accretion rates ob- 
tained from the simulations. However, since their fitting 
formula is not based on cosmology-independent quanti- 
ties, it is only valid for the particular cosmology adopted 
in the Millennium simulations (see Section 15.11 for more 
details) . 

4.3. Redshift distribution of subhalo accretion 

In Fig. [7] we show our model predictions for the 
distribution of the accretion redshift for subhalos with 
m a /M = 0.1 (solid lines), 0.03 (dotted lines) , 0.01 
(dashed lines), 0.003 (long dashed lines) and 0.001 (dot- 
dashed lines), respectively. The results are obtained us- 
ing Eq. [S] by making some simple coordinate transfor- 
mations. Results shown in different panels are for host 
halos of different masses, as indicated. Open circles and 
filled triangles indicate the results obtained from the 
300 h~ 1 Mpc and 100 /i _1 Mpc simulations boxes, respec- 
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./V-body simulations with the same cosmology (open circles). For comparison, results obtained from the 100 h~ 1 Mpc box simulations are 
also shown (as filled triangles) for cases where statistics are sufficiently good. 



tively, where the error-bars have been obtained using 200 
bootstrap resamples. The various lines show the predic- 
tions based on Model III, and overall match the simu- 
lation results remarkably well. Note that the accretion 
rate depends strongly on the mass of the host halo. For 
the same mass ratio, subhalos in more massive hosts are 
accreted later, reflecting the hierarchical nature of struc- 
ture formation in the ACDM cosmology. 

4.4. Un-evolved subhalo mass functions 

Finally, let us look at the un-evolved subhalo mass 
functions. By integrating Eq. ([3]) over a given redshift 
range, we can obtain the un-evolved mass function of 
the subhalos accreted in that redshift range. In Fig. |8] 
we show the un-evolved mass functions of subhalos ac- 
creted in the redshift ranges [0, 1], [1, 2], [2, 3], [3, 4] 
and [4, 5], respectively. Results are shown for host ha- 
los of different masses, as indicated in each panel. Here 



again, symbols indicate the results from our simulation 
boxes, while lines show the predictions of Model III. 
Clearly, our model is in excellent agreement with the 
simulation results at all redshifts and for all host masses. 
Upon close inspection, it is clear that the un-evolved sub- 
halo mass function for a given redshift range depends on 
host halo mass, especially at high redshift: in terms of 
the scaled mass, m a /Mo, the subhalo mass function at 
high z is significantly higher for lower-mass host halos. 
Moreover, the normalization of the un-evolved subhalo 
mass function at a given redshift for halos of different 
masses seem to be roughly proportional to the assem- 
bly history of the host halos shown in Fig. [TJ To test 
this, we show in Fig. [9] the un-evolved subhalo mass 
functions for different host halos at the time when the 
host halos have assembled a fixed fraction of their fi- 
nal masses, i.e. for subhalos accreted in a given range 
of \og[A4 a /M ] range. Results are shown for five dif- 





log[m a /M ] log[m a /M ] 

Fig. 8. — The un-evolved mass function of subhalos accreted within the redshift ranges [0, 1] (solid lines), [1, 2] (dotted), [2, 3] (dashed), 
[3, 4] (long dashed) and [4, 5] (dot-dashed), for host halos of different masses, as indicated in each panel. Here model predictions assuming 
a ACDM cosmology are compared with the the results obtained from the 300 /i — x Mpc box ./V-body simulations of the same cosmology 
(open circles). For comparison, results obtained from the 100 h~ *Mpc box simulation are also shown (as filled triangles) for cases where 
the statistics are sufficiently good. The vertical lines in each panel correspond to the mass of 100 particles in the two simulations, as in 
Fig. HJ. 



ferent ranges of log[M a /M ]: [—0.5,0.0] (solid lines), 
[-1.0,-0.5] (dotted), [-1.5,-1.0] (dashed), [-2.0,-1.5] 
(long dashed), and [—2.5,-2.0] (dot-dashed). For each 
range of log[A^ a /A/ ] the subhalo mass functions for 
four host halo masses, M Q ~ 10 145 , 10 14 , 10 13 , and 
1O 12 '°/i _1 M0, are plotted with the same line style. In- 
terestingly, the mass function of subhalos accreted in a 
given range of \og[M a /M } is almost independent of the 
host halo mass, demonstrating that the amplitude of of 
the un-evolved subhalo mass function at a given redshift 
is directly related to the mass assembly rat e of the host 
halo at that redshift (see also lGiocolill2010l ). 

Integrating Eq. ([3]) over the full redshift range (here 
for practice we adopt redshift range z = — 10) yields 
the total un-evolved subhalo mass function. The results 
are plotted in the left-hand panel of Fig. [TO] for host ha- 



los of different masses. Interestingly, although host ha- 
los of different masses accrete their subhalos at different 
redshifts, their total un-evolved subhalo mass functions 
are extremely similar, both in the simulations (symbols) 
and in the model (lines). Note, though, that this uni- 
versality of the un-evolved subhalo ma ss function, first 
hinted at in Ivan den Bosch et al.l (|2005[ ). is only approxi- 
mate. This is due to the fact that the un-evolved subhalo 
mass function depends on the shape of the perturbation 
power spectrum around the mass scale in question. For 
the ACDM cosmology considered here, the slope of the 
power spectrum only changes by a small amount over 
the mass range 10 12 /i^Mq < M < 10 15 h. -1 M . In- 
deed, close inspection of the left-hand panel of Fig. [TO] 
reveals small differences between the curves for different 
host masses. 
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Fig. 9. — The un-evolved mass function of subhalos accreted 
within different log [Ala /Mo] ranges [-0.5, 0.0] (solid lines), [-1.0, 
-0.5] (dotted), [-1.5, -1.0] (dashed), [-2.0, -1.5] (long dashed) and 
[-2.5, -2.0] (dot-dashed). For each range of Iog[jVf a /Mo], results 
are shown with the same line style for host halos of four different 
masses: M ~ 10 145 , 10 140 , 10 130 , and 10 12 - /l _1 M Q . Note 
that shown in this way the results for different Mo are almost 
indistinguishable. 

To further illustrate the dependence of the un-evolved 
subhalo mass function on the shape of the power spec- 
trum, the right-hand panel of Fig. [TU] shows the model 
predictions for the total un-evolved subhalo mass func- 
tions for scale-free models with power spectra with spec- 
tral indices n — (dashed line), n = — 1 (long dashed 
line) and n = — 2 (dot-dashed line), respectively, together 
with that for 10 12 /i _1 M Q host halos in our ACDM cos- 
mology (dotted line). The prediction for the n = — 2 
model is very close to that of the ACDM model, reflecting 
the fact that the effective spectral index of the ACDM 
power spectrum is close to —2 over the mass scales in 
question. The subhalo mass functions for the n — and 
n = — 1 models are much shallower. 

Finally, for comparison, the solid line in the left-hand 
panel of Fig. [TU1 show s the fitting formula obtained by 
iGiocoli et al.l (|2008a[ ) from A^-body simulations of the 
ACDM model. This fitting formula agrees well with our 
model prediction for the ACDM model, except at the 
high mass end where it is slightly lower than our model 
prediction. Compared to the N-body simulation results, 
our model prediction agrees with the data slightly better, 
except perhaps at the very high massive end. 

5. DISCUSSION 

In this paper we have developed an analytical model 
for the mass function of CDM subhalos at their time 
of accretion and for the distribution of their accretion 
times. This model can be used to predict the un-evolved 
subhalo mass function, the mass function of subhalos ac- 
creted at a given time, the accretion-time distribution 
of subhalos of a given initial mass, and the frequency of 
major mergers as a function of time. We have tested our 
model against results obtained from high-resolution N- 
body simulations, and found that the model predictions 
match the simulation results extremely well. In this sec- 



tion, we first discuss the universality of our model based 
on the ingredients used in the construction of the model. 
We then briefly describe two possible applications of our 
model. 

5.1. The Universality of the Model 

The first ingredient of our model is the adopted form 
for F(s a , 8 a \So, Sq] M a ). As described in Section [37T1 
all of our three models, I, II and III, are based on the 
extended PS theory, which has been shown to work 
reasonably well for various hierar chical models (e.g. 
lLacev fc CoIiH99l iCole et al.ll2008l ). The modifications 
made in Models II and III have not been tested for other 
cosmological models. However, the fact that these mod- 
els work well for the ACDM model over a large redshift 
range, during which the cosmological parameters change 
by fairly large amounts, suggests that the modifications 
should also work for other cosmological models with sim- 
ilar power spectra. Unfortunately, the effective power in- 
dices covered by the ACDM spectrum are limited, and it 
is presently unclear whether the modifications will work 
equally well for models with vastly different power spec- 
tra. It will be interesting to test the validity of our model 
using scale-free models that cover a large range in spec- 
tral indices. 

The second ingredient of our model is the form of 
P(M„.\St) ,St)) described in Section [3T2l As tested exten- 



sively in IZhao et al.l (|2009( ). their model for the median 
accretion history is universal and works not only for re- 
alistic CDM models but also for scale-free models with 
different spectral indices. Since our tests cover a range of 
rcdshifts over which the cosmological parameters change 
by a large amount, we expect the log- normal form of 
P(M a \So, So), and the mass dependence of its dispersion 
(Eq. [21]), to also apply to other models that have a 
power-spectrum that is not too different from that of the 
ACDM cosmology considered here. Here again, it will 
be interesting to test the validity of our model for other 
power spectra using scale-free models. 

Similar arguments also apply to the correction given 
by Eg. (|22| . This ' correction' is empirically obtained 



by IZhao et al.l (pOOl ) from various CDM and scale-free 
models. As discussed in £)3.2l using the model with or 
without this 'correction' corresponds to two different def- 
initions of the un-evolved subhalo population, and only 
has a non-negligible impact when it comes to subhalos 
accreted at low redshifts in low-mass host halos. 

To summarize, we believe that our model holds for 
any variant of the ACDM model. In particular, it should 
work accurately for any cosmological model whose pa- 
rameters are in agreement with current observational 
constraints from the cosmic microwave background, su- 
pernovae la, and galaxy redshift surveys. 

5.2. Applications 

Our model has a number of applications. Here we focus 
on two of them: (i) the evolution of subhalos as they 
orbit within their hosts, and (ii) the construction of a 
self-consistent model to link galaxies and dark matter 
halos across cosmic time. 

The model described here applies to subhalos at ac- 
cretion, i.e. the un-evolved population of subhalos. 
However, as a small halo merges with and orbits in 
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Fig. 10. — Left panel: the un-evolved subhalo mass distribution as a function of m a /Mo for host halos of different final masses at z = 0. 
Here model predictions (lines) are compared with the results obtained from both the 300 h~ x Mpc box (open circles) and 100 h~ x Mpc box 
(filled triangl e s) iV-b ody simulations of the same cosmology. For comparison, we also include in the figure the fitting formula obtained by 
Giocoli ct al. (2008a) as the solid line. Right panel: The predicted un-evolved subhalo mass functions for scale-free models with n = 0, —1 
and —2. For comparison, we also include in the figure the model prediction for the ACDM model considered in this paper. 



a larger halo, it is subjected to tidal forces from the 
host, causing it to lose mass, and to dynamical friction, 
which causes it to lose energy and angular momentum 
to the dark matter particles of the host. Thus, dy- 
namical evolution after accretion can change the prop- 
erties of the subhalo population. A great deal of 
work has been done to understand the properties of 
the 'evolved' subhalo population, using e ither numer- 
ical s imulations (e.g. [Klvpin et al.l 119991: iMoore et all 
19991 : iSpringel et aLN2001bt iDiemand et al.l 12004 120071 



Gao et all 12004 Kang et all 120051: iSpringel et al] 12008 
Wetzel fc W hite 20 13), or merger trees constructed from 



the e x tende d PS formalism (e.g. [Taylor fc Babul I2001L 
2001 12001 Ivan den Bosch et all 120051: IZentner et all 
20051 IGiocoli et alll2008al: iGan et al.ll20ldl 

After accretion, whether a subhalo can survive as 
a self-bound entity depends on its mass (relative to 
that of its host), its density profile (which is related to 
its accretion redshift), and its orbit. The model pre- 
sented here provides information about the first two 
parts. Thus, combined with infor mation about the ini- 
tial o r bits of accreted subha l os (e.g . Ivan den Bosch et"all 
19991: iKhochfar fc Burkertl 12001 



Ludlow et all 120091: 



Valluri et al.ll2010l : lWetzelll201lD and about how dynam- 



ical friction and tidal stripping operate on su bhalos (e.g. 
IJiang et al.ll200l IBovlan-Kolchin et al.ll2008D . our model 
can be used to construct models for the evolved subhalo 
population, such as their mass function, their distribu- 
tion in host halos, and their correlation with host halo 
properties. We will come back to such modeling in a 
future paper. 

The last decade has seen much effort in using 
halo occupation statistics to describe the galaxy-dark 
matter connection in an attempt to understand the 
galaxy luminosity function, galaxy clustering, galaxy- 
galaxy l ensing, and the kinemati c s of satellite galax- 
ies (e.g . IMo. Mao fe White! 119991: IBerlind fe Weinberg! 

20031: Ivan den Bosch et all 1 2003 1 : 
Tin ker et all 120051: iCooTavl [2005, 



2 0021; rfang et all 
Zheng et al.M 2005: 



2006 


; 

: 


Coorav & Ouchil 120061: IVale & Ostriker 




2004, 


2006 


van den Bosch et all 120071; Yang ct al. 


2009; 


Cacciato et all 120091: iMore et all l2009t iMoster et al. 



a direct link between the halo occupation statistics 
of galaxies and of dark matter subhalos, using abun- 
dance matchi ng techniques to link satellite galaxies to 
subh al os (e.g. IKravtsov et al.l 12004 IConrov et all 1 20061. 



20071: iConroy fe Wechslerl 120091: iBehroozi et al.l I20K 
Wang fe JingfeOlOHWetzel "fe Whitd l2010HNeis tein et all 
201 lUAvila- Reese fc Firmanil l20l!l) . These investiga- 
tions typically rely on the un-evolved subhalo mass func- 
tion, but do not account for the possibility that differ- 
ent subhalos are accreted at different times and that 
the halo mass - galaxy mass relation may be redshift- 
dependent. In order to construct a self-consistent model 
based on abundance matching, one needs not only the 
(un-evolved) subhalo mass function, but also the distri- 
bution of accretion times. This is exac tly what our mode l 
can provide. In a forthcoming paper (|Yang et alll201ll) . 
we will construct such a self-consistent model to char- 
acterize how the galaxy-dark matter connection evolves 
over time. 
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